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1  Summary  of  Developments 

This  report  documents  activities  completed  under  AFOSR  Grant  F49620-02-1-0068 
from  January  2002- June  2005.  Although  activities  under  this  grant  were  completed 
in  December  2004,  work  towards  its  original  objectives  was  only  recently  completed 
on  June  30th  2005.  In  particular,  some  of  the  research  directions  that  were  initially 
pursued  could  only  be  completed  on  that  date.  Three  graduate  students  were  each 
partially  supported  on  this  grant,  one  of  whom  (Dr.  Debraj  Ghosh)  successfully 
defended  his  doctoral  thesis,  at  the  Johns  Hopkins  University,  on  June  30  2005,  while 
the  other  two  students  (Mr.  Sonjoy  Das  and  Mr.  Alireza  Doostan)  have  already 
completed  significant  work  towards  their  dissertations.  Their  work  on  this  AFOSR 
portion  of  their  work  will  culminate  in  a  MS  thesis  for  each  of  them. 

The  activities  under  this  grant  were  in  the  general  area  of  Verification  and  Val¬ 
idation  using  the  Polynomial  Chaos  formalism  developed  by  the  PI  over  the  course 
of  the  past  15  years.  In  particular,  three  significant  questions  were  formulated  and 
addressed  in  the  course  of  this  research: 

1.  How  to  develop  polynomial  chaos  representations  of  physical  parameters  from 
experimental  measurements  of  these  parameters  ?  This  provides  the  link  be¬ 
tween  experimental  evidence  and  the  stochastic  predictive  process.  It  has 
enabled  us  to  eliminate  severe  assumptions  associated  with  probabilistic  ap¬ 
proaches  in  general,  namely  having  to  invoke  various  assumptions  as  to  the 
probabilistic  measure  of  the  data  (gaussian,  lognormal,  various  Askey  scheme 
measures,  ...) 

2.  How  to  propagate  the  uncertainty  in  these  parameters  (as  reflected  in  their 
Polynomial  Chaos  representations)  into  the  dynamical  behavior  of  the  physi¬ 
cal  system.  Specifically,  a  novel  approach  is  developed  to  the  random  eigen¬ 
value  problem  using  Polynomial  Chaos  representations  coupled  with  Stochastic 
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Galerkin  Projections.  The  new  approach  provides  both  efficient  algorithms  for 
characterizing  the  solution  to  the  random  eigenvalue  problem  as  well  as  new 
insight  as  to  the  significance  and  usefulness  of  this  solution. 

3.  How  important  is  data  refinement:  what  is  the  significance,  on  the  predictive 
value  of  a  computational  model,  of  collecting  additional  experimental  measure¬ 
ments  as  opposed  to  performing  further  numerical  refinement  (using,  for  exam¬ 
ple,  mesh  refinement).  This  ’’data  refinement”  issue  is  significant  if  predictive 
tools  are  to  become  surrogates  to  physical  experiments. 

The  remainder  of  this  report  will  review  the  highlights  of  each  of  the  above  three 
topics.  Details  of  the  analysis  are  contained  in  seven  papers  which  are  various  stages 
of  the  publication  process:  one  is  in  press  for  the  AIAA  Journal  [24],  one  has  been 
accepted  for  publication  in  the  International  Journal  of  Numerical  Methods  in  Engi¬ 
neering  [7],  one  is  under  review  for  the  IJNME  [18,  23],  one  is  under  review  for  the 
Journal  of  Computational  Physics  [17],  and  two  are  still  in  preparation  [8,  15]. 

2  Characterization  of  Chaos  Representations  from 
Measurements 

In  a  series  of  previous  publications,  the  PI  presented  a  mathematical  framework  for 
the  characterization  and  propagation  of  uncertainty  in  physical  systems  [31,  16,  13, 
12,  14,  27,  19,  22,  21,  20,  29,  28,  6,  5,  24],  That  work  is  based  on  the  adaptation 
of  multiple  Wiener  integral  representations  [33,  4]  to  finite-dimensional  spaces  and 
their  implementation  into  a  weighted  residual  scheme  for  the  stochastic  characteri¬ 
zation  of  the  solution  of  stochastic  partial  differential  equations.  The  restriction  of 
the  representations  to  finite-dimensional  uncertainty  (i.e.  stochastic  processes  char¬ 
acterized  by  a  finite-dimensional  random  vector)  permits  the  generalization  of  the 
Wiener  constructions  which  had  used  polynomials  orthogonal  with  respect  to  Gaus¬ 
sian  measure.  Extensions  of  that  work  to  the  Askey  scheme  were  recently  developed 
[35,  34]  together  with  extensions  using  non-orthogonal  representations  [1,  2]  and  rep¬ 
resentations  in  terms  of  wavelets  [25,  26].  While  the  above  extensions  are  limited 
to  representations  in  terms  of  independent  random  variables,  extensions  using  finite¬ 
dimensional  dependent  random  vectors  (such  as  appearing,  for  instance,  in  the  finite¬ 
dimensional  Karhunen-Loeve  representation  of  an  arbitrary  stochastic  process)  have 
also  been  completed  [30].  These  so-called  polynomial  chaos  expansions,  coupled  with 
stochastic  projection  mechanisms,  provide  a  general  method  for  characterizing  the 
solution  of  problems  of  mathematical  physics  whose  parameters  have  been  described 
as  stochastic  process.  An  analysis  of  the  error  associated  with  this  method  has  al¬ 
ready  been  developed  [3,  2].  The  mathematical  foundation  of  these  developments, 
which  permits  their  rigorous  mathematical  analysis,  lies  in  functional  analysis  and  in 
the  observation  that  second-order  random  vectors  (i.e.  loosely  speaking,  those  with 
finite  second  order  statistics)  form  a  Hilbert  space. 
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Implicit  in  all  these  developments  is  the  assumption  that  the  parameters  in  the 
governing  equations  have  been  accurately  characterized  as  stochastic  processes.  In 
this  case,  the  formalism  described  above,  based  on  the  Chaos  expansions,  can  be 
viewed  as  an  efficient  procedure  to  propagate  the  probabilistic  measures  from  the 
system  parameters  to  the  solution,  and  as  such  appears  to  be  an  alternative  to  ap¬ 
proaches  based  on  Monte  Carlo  sampling  or  perturbation  expansions  [36].  As  noted 
above,  it  is  often  the  case  that  system  parameters  are  not  known  with  enough  res¬ 
olution  to  permit  their  accurate  characterization  as  stochastic  processes.  In  such 
cases,  the  approach  based  on  Chaos  developments  provides  a  unique  perspective  on 
the  problem  as  the  sensitivity  of  the  system  parameters  to  additional  information  can 
be  cast  as  perturbing  their  coordinates  with  respect  to  the  polynomial  chaos,  which 
in  turn  is  readily  described  as  perturbations  in  Hilbert  space.  The  impact  of  these 
perturbations  on  the  chaos  coordinates  of  the  solution  to  the  governing  equations  can 
then  be  viewed  as  quantifying  the  impact  of  refining  the  probabilistic  measure  of  the 
data  on  the  predictive  capability  of  the  mathematical  model.  Ours  efforts  in  this 
direction  have  relied  on  the  maximum-likelihood  arguments  to  compute  estimates  of 
the  Chaos  coordinates  of  parameters  from  associated  statistical  samples  [7].  The  same 
problem  has  also  been  addressed  using  the  Bayesian  framework  for  parameter  estima¬ 
tion.  The  benefit  of  this  approach  lies  in  its  ability  to  characterize  the  statistics  of  the 
estimates,  thus  enabling  the  determination  of  their  accuracy  and  their  sensitivity  to 
further  data.  In  addition,  a  propagation  of  the  error  associated  with  these  estimates 
to  the  predictions  of  the  stochastic  governing  equations  has  been  developed.  This 
clearly  provides  an  essential  ingredient  for  any  model  validation  process. 

The  procedure  described  in  this  section  is  described  in  various  recent  publications 
by  PI  and  co-workers  [9,  17,  7] 

3  Data  Refinement 

By  representing  the  coefficients  in  the  Chaos  representations  as  random  variables 
themselves,  it  becomes  possible  to  determine  the  influence  of  the  uncertainty  in  their 
value  on  the  predicted  solution. 

Specifically,  as  experimental  data  is  obtained,  estimates  are  updated  for  any  statis¬ 
tics  of  the  stochastic  process  being  represented.  In  particular,  the  chaos  coefficients 
are  a  particular  kind  of  statistics.  When  a  stochastic  representation  is  available  for 
these  coefficients,  reflecting  the  uncertainty  in  their  value,  the  impact  of  this  uncer¬ 
tainty  can  be  propagated  to  the  solution,  and  the  worth  of  additional  information 
can  be  ascertained.  It  is  assumed  that  the  uncertainty  in  the  chaos  coefficients  is  in¬ 
dependent  of  the  uncertainty  in  the  mechanistic  parameters  themselves.  This  results 
in  an  increase  in  the  dimension  of  the  chaos  in  which  the  solution  is  described. 

The  problem  then  becomes  a  standard  problem  in  stochastic  computational  me¬ 
chanics,  and  can  be  very  efficiently  solved  using  standard  procedures  established  by 
the  PI  for  Stochastic  Galerkin  Projections.  It  should  be  noted  here  that  the  so-called 
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’’generalized  chaos”  expansions  in  terms  of  non-gaussian  measures  cannot  be  applied 
to  this  task  without  significant  assumptions,  specifically  in  terms  of  completeness  in 
infinite  dimensional  spaces  (in  this  case,  the  tensorized  construction  of  bases  using 
one-dimensional  spaces  as  building  blocks  only  works  for  the  Gaussian  and  Poisson 
measures  [32,  10,  11]). 

The  procedure  described  in  this  section  is  described  in  an  invited  paper  to  the 
Journal  of  Computational  Physics,  and  has  been  presented  in  conference  proceedings 
[9,  8,  17] 

4  Random  Eigenvalue  Problem:  Polynomial  Chaos 
Representation 

The  polynomial  chaos  decomposition  involves  representing  a  stochastic  process  with 
respect  to  a  Hilbertian  basis  consisting  of  the  multidimensional  Hermite  polynomi¬ 
als  of  orthogonal  gaussian  variables.  For  the  random  eigenvalue  problems  where  the 
stiffness  K  is  a  function  of  a  set  of  random  variables  £,  the  polynomial  chaos  repre¬ 
sentation  of  the  eigenvectors  and  eigenvalue  has  also  been  used.  Thus  starting  with 

K(£)4>  =  XM<j>  (1) 

The  eigenvalues  and  the  eigenvectors  are  then  expressed  as 

p  p 

(t>  =  fafa,  a  =  fafa 

i=0  i=0 

where  fa  are  zero-mean,  orthogonal  polynomials  in  £  such  that, 

fa  =  1;  (fa)  =0  i  >  0;  (fa fa)  =  $ij(fal)  i,j  >  0. 

fa  and  A;,  the  coefficients  in  the  expansion,  are  calculated  using  their  generalized 
Fourier  coefficient  expressions, 

jt-  _  (4>fa)  ,  i^fa) 

^  m  ’  * "  m  ' 

The  denominator  in  the  above  expressions  can  be  exactly  evaluated,  whereas  the 
numerator  must  be  estimated,  for  example,  by  Monte  Carlo  sampling.  In  this  case, 
each  realization  of  the  set  of  random  variables  £  is  associated  with  a  specific  realization 
of  the  stochastic  system  and  the  corresponding  eigenproblem  is  solved  for  A  and 
4>.  The  realization  of  fa  is  simultaneously  computed  as  a  polynomial  form  in  the 
set  {£j}.  The  coefficients  fa  and  A,  are  then  estimated  by  interpreting  the  inner 
product  as  a  statistical  average  estimated  using  the  realizations.  Once  the  chaos 
coefficients  have  been  computed,  realizations  of  the  eigenvalues  or  eigenvectors  can 
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be  readily  synthesized  and  used  to  estimate  higher  order  statistics  of  the  probability 
density  functions.  In  particular,  the  first  and  second  order  statistical  moments  of  the 
eigenvalues  are  obtained  as, 

p 

A  =  (A)  =  A0,  Var(A)  =  £>?)A ? 

i=  1 

Similarly,  for  the  eigenvectors, 

<i>  =(<}>)  =  4>q  ■ 

It  can  be  noted  that  the  polynomial  chaos  approach  does  not  require  explicit  knowl¬ 
edge  of  the  derivatives  of  the  stiffness  and  mass  matrices  with  respect  to  the  random 
parameters. 

The  procedure  described  in  this  section  has  been  presented  in  journal  publications 
[23,  18]  and  various  conference  proceedings. 

4.1  Intrusive,  Galerkin-Based  Approach 

An  alternative  way  to  compute  the  coefficients  in  the  chaos  representation  of  the 
eigensolution,  which  does  not  require  Monte  Carlo  sampling,  is  to  pose  the  problem 
as  an  optimization  problem.  In  particular,  the  chaos  coefficients  of  eigenvalues  and 
eigenvectors  are  sought  such  that  the  error  between  the  left  and  right  sides  of  the 
eigen-equation  is  minimized.  This  is  described  next. 

Consider  the  random  eigenvalue  problem  given  by  the  equation 

K(0<t>  =  \<t>  ■  (2) 

Let  also  expand  K  in  polynomial  chaos  expansion, 

L— 1 

K  =  Yl^Ki-  (3) 

i=0 

Note  that  an  L-term  expansion  is  used  for  the  matrix  K  representing  the  data,  and 
a  P-term  expansion  is  used  for  the  unknown  eigenvectors  and  the  eigenvalues.  Now 
onward,  the  subscript  l  is  removed  from  A  snd  4>  for  convenience,  implicitly  assuming 
that  the  rest  of  the  formulation  is  done  for  a  single  mode.  Let  us  define  Substituting 
Chaos  expansions  into  the  eigen-equation  results  in 

L-l  P-i  P—l P—l 

Y  Y  ^iKi<pb)  =  Y  S  a(,vo)  ■  (4) 

i=0  j=0  i=0  j— 0 

The  above  equality  is  next  interpreted  to  hold  in  the  weak  sense.  Thus  multiplying 
the  equation  by  ipk  and  taking  the  inner  product  in  the  Hilbert  space  of  second  order 
random  variables  results  in, 
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A:  =  0, . . . ,  P  —  1  ,  (5) 


L-1P-1  P-1  P-1 

EE  EE 

i= 0  j=0  i= 0  j= 0 

that  can  be  rewritten  as, 


K$>  =  A$  .  (6) 

Here  $  is  nP  dimensional  vector,  constructed  by  arranging  P  number  of  n  dimensional 
vectors  <f>^  in  a  single  column.  K  and  A  are  the  matrices  of  dimension  ( nP  x  nP) , 
described  below,  set  of  vectors  and  scalars  describe  the  behavior  of  a  single  physical 
mode  only. 

A  Newton  method  is  used  to  solve  the  associted  optimization  problem. 

It  should  be  noted  that  in  the  Monte  Carlo  simulation,  and  in  the  existing  sim¬ 
ulation  based  chaos  expansion  method,  the  eigenvalues  are  ordered  in  ascending  or 
descending  order  according  to  their  numerical  values  in  each  realization.  For  repeated 
eigenvalue  cases  and  some  closely  spaced  eigenvalue  cases  this  ordering  strategy  results 
in  modal  switching  whereby  the  relative  ordering  changes  for  two  modes  associated 
with  two  distinct  behavior.  The  closeness  of  the  modes  are  measured  with  respect 
to  the  strength  of  perturbation.  This  Galerkin-based  technique  deals  with  statis¬ 
tical  description  of  both  the  eigenvalues  and  eigenvectors  simultaneously,  resulting 
in  eigenpair  ordering  rather  than  the  eigenvalue  ordering.  These  two  ordering  tech¬ 
niques  differ  in  the  above  mentioned  closely  spaced  or  repeated  eigenvalues  cases.  In 
such  cases,  studying  the  behavior  of  the  modal  subspace  could  be  more  useful  rather 
than  considering  the  individual  modes.  However,  this  paper  considers  only  the  widely 
spaced  modes  where  such  situations  does  not  arise. 

The  non-intrusive  approach  is  described  in  detail  in  [24]. 

4.2  Application  to  large  scale  system 

Implementing  the  present  concepts  in  the  context  of  large-scale  systems  (millions  of 
degrees  of  freedom)  presents  a  number  of  serious  challenges.  A  number  of  issues  have 
been  encountered  and  resolved. 

1.  The  size  of  the  jacobian  arising  in  the  Newton’s  method  is  very  large,  thus  using 
the  direct  linear  solvers  become  prohibitive,  an  iterative  solver  is  needed.  For 
a  small  system,  it  is  observed  that  the  MINRES  algorithm  is  a  very  efficient 
choice. 

2.  Mat-vec  multiplication:  Due  to  the  large  size  of  the  jacobian,  it  can  not  be 
stored,  thus  all  the  rows  are  recomputed  whenever  a  matrix- vector  product  is 
sought.  This  results  in  heavy  computational  burden.  More  efficient  methods 
for  the  mat-vec  product  are  needed,  this  is  a  crucial  hinge. 
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4.3  Repeated  and  Closely  Spaced  Modes 

If  the  eigenvalues  are  clustered,  then  the  corresponding  eigenvalues  and  the  eigen¬ 
vectors  are  highly  sensitive  to  the  system  parameters.  In  the  chaos  expansion  of 
such  eigenvalues  and  the  eigenvectors,  two  problems  are  observed.  First,  neither  the 
non-intrusive  nor  the  intrusive  methods  can  adequately  capture  the  variability  of  the 
eigensolution.  Second,  the  results  obtained  using  these  two  methods  are  quite  differ¬ 
ent.  But  both  of  them  can  capture  the  stable  invariant  subspace  formed  by  the  set  of 
the  eigenvectors  corresponding  to  the  cluster  of  the  eigenvalues.  An  analysis  of  the 
problem  can  be  done  using  a  simple  mass-spring  system. 

Recently  an  enriched  version  of  the  polynomial  chaos  expansion  has  been  proposed. 
The  idea  is  to  simply  augment  the  set  of  basis  functions  with  some  special  functions 
that  anticipate  particular  behaviors  in  the  response  of  the  system.  These  functions 
will  typically  consist  of  a  step  function,  an  absolute-value  function,  or  a  rational 
polynomial.  For  a  system  with  one  random  parameter,  it  is  observed  that  for  the 
expansion  of  the  eigenvalue,  adding  an  absolute- value  function  and  for  the  expansion 
of  the  eigenvector,  adding  a  step  function  improves  the  result  significantly. 

Once  the  Chaos  representation  has  been  computed,  a  corresponding  large  statis¬ 
tical  sample  can  be  easily  generated,  from  which  estimates  can  be  developed  for  the 
probability  density  function  of  the  response  quantity  of  interest. 

Results  from  this  part  of  the  work  have  been  presented  at  various  conferences  and 
submitted  for  publication  in  a  refereed  journal  [23]. 

5  On-Going  Work 

Dr.  Ghosh  is  spending  the  next  two  months  completing  the  implementation  of  the 
stochastic  eigen-analysis  to  large  scale  and  arbitrary  systems.  He  is  supported  on  this 
work  by  a  discretionary  account  available  to  the  PI  at  USC. 

Mr.  Doostan  will  be  finishing  his  Doctoral  thesis  within  the  next  15  months.  He 
is  working  on  integrating  the  errors  from  data  limitations,  numerical  truncation  and 
statistical  sample  limitations  in  order  to  formulate  a  comprehensive  error  budget  that 
can  be  used  in  defining  a  rigorous,  quantifiable  and  constructive  validation  process. 
Mr.  Doostan’s  research  is  supported  by  a  recent  NSF  grant. 

Mr.  Das  will  be  also  be  finishing  his  doctoral  thesis  within  the  next  15  months. 
He  is  working  on  alternative  procedures  for  polynomial  chaos  characterization. 
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